function plot_sol(sol, nx, ny)
an = zeros((nx+1),(ny+1));
it = 1;
for j = 1:ny+1
    for i = 1:nx+1
        an(i,j) = sol(it);
        it = it +1;
    end
end
it
size(sol)
sol
figure(1);
imagesc(an);
colorbar;
%colormap(grey);

figure(2);
surf(an);
figure(3);
contour(an);